
import numpy as np
import math
from scipy import linalg


#test
scalex = np.array([[31.0,-13,0,0,0,-10,0,0,0],[-13,35,-9,0,-11,0,0,0,0],[0,-9,31,-10,0,0,0,0,0],[0,0,-10,79,-30,0,0,0,-9],
                [0,0,0,-30,57,-7,0,-5,0],[0,0,0,0,-7,47,-30,0,0,],[0,0,0,0,0,-30,41,0,0],[0,0,0,0,-5,0,0,27,-2],
                [0,0,0,-9,0,0,0,-2,29]])
scaley = np.array([[-15.0],[27],[-23],[0],[-20],[12],[-7],[7],[10]])
testx = np.array([[1,2,3,4],[2,9,12,15],[3,26,41,49],[5,40,107,135]])
testy = np.array([1,2,3,4])
scalex2 = np.array([[30,33,-43,-11,-38,-29,37,28,23],[-480,-523,644,128,621,480,-618,-489,-329],
                    [60,266,-1862,-1991,464,546,-968,-1567,1652],[540,624,-782,290,-893,123,567,5,-122],
                    [-450,-675,2245,2326,-1512,1230,-822,129,-189],[-300,-120,-1114,-1295,1946,302,-376,-1540,-609],
                    [1080,998,508,2460,-1628,-1358,2896,2828,-2002],[-1080,-1408,3340,2267,21,-1202,866,-2690,-1351],
                    [-300,-435,1594,1685,340,2279,-27,2917,-2336]])
scaley2 = np.array([188,-3145,-4994,680,7845,1876,9712,-11599,10127])
scalex = [[45,80,34,56,98,56,24,74,52,92],[66,74,74,34,72,46,77,28,84,27],[47,33,23,34,31,72,85,77,74,63],
                  [40,22,78,63,97,40,28,96,53,91],[31,24,66,59,96,96,36,77,34,23],[44,100,67,34,38,27,84,54,35,100],
                  [25,97,99,95,45,43,70,97,46,97],[35,72,55,28,52,81,92,39,74,55],[34,61,99,51,38,20,56,66,69,68],
                  [72,27,50,48,51,78,61,56,96,37]]
scaley = [62,58,64,96,84,90,74,34,40,25]

#换元  表示传入的x矩阵,y矩阵,第m行,总的行数n
#print(scalex2[0,:])
#----------------参数------------------
n = 20
import copy 
np.random.seed()

rd = np.random.RandomState(888)
radomx = rd.uniform(1,100,(n,n))
radomx2 =copy.deepcopy (radomx)
#print(radomx)
#print(radomx2)
radomy = rd.uniform(100,1000,(n,1))
radomy2 = copy.deepcopy(radomy)
#print(radomy)
def swap(scalex,scaley,m,n):

    temp = scalex[m][m]
    maxcowls = m
    for i in range(m,n):
        if(temp<abs(scalex[i][m]) ):
            temp = abs(scalex[i][m])
            maxcowls = i
          
    
    tempscalex = np.copy(scalex[m])
    tempscaley = np.copy(scaley[m])
    scalex[m] = scalex[maxcowls,:]
    scalex[maxcowls,:] = tempscalex
    scaley[m] = scaley[maxcowls]
    scaley[maxcowls] = tempscaley



#列主元素消元法
def gaussinlist(scalex,scaley):
    m= len(scalex)
    n = m
    testx = copy.deepcopy(scalex)
    if(m>n):
        print('请输入正确的矩阵')
        return False
    else:
        #初始化参数矩阵
        l = np.zeros((n,n))
        for i in  range(0,n):
            if(scalex[i][i]==0):
                print('主元素为0')
        for k in range(n-1):
            swap(scalex,scaley,k,n)
            for i in range(k+1,n):
                #计算每次需要乘的数
                l[i][k] = scalex[i][k] / scalex[k][k]
                for j in range(n):
                    try:
                        scalex[i][j] =round((scalex[i][j] - (l[i][k] * scalex[k][j])),2)
                    except:
                        print(scalex[k][i],"i = ",i)
                        print(scalex[j][i])
                        print(l[k][j])
                        exit()
                scaley[k] = round((scaley[k] - l[k][j] * scaley[j]),2)
                print(scaley)
    #形成三角矩阵,求解方程
        print(scalex)
        print(scaley)
        x = np.zeros(n)
        y = np.zeros(n)
        x[n-1] = scaley[n-1] / scalex[n-1][n-1]
        for i in range(n-2,-1,-1):
            for j in range(i+1,n):
                scaley[i] = scaley[i]- (x[j] * scalex[i][j])
            x[i] = scaley[i] / scalex[i][i]
        for i in range(0,n):
            print("x" + str(i+1) + "=" ,x[i])
        for i in range(0,n):
            for j in range(0,n):
                y[i] += testx[i][j]*x[j]
       # print(y)        
def SolvingLinearEquations(AugmentedMatrix):
    col = AugmentedMatrix.shape[1] #增广矩阵列数
    #消元
    for i in range(col- 2):
        current_column  = AugmentedMatrix[i:,i] 
        max_index = np.argmax(current_column) + i #寻找最大元
        if(AugmentedMatrix[max_index,i] == 0):
            print("无唯一解")
            return
        AugmentedMatrix[[i,max_index],:] = AugmentedMatrix[[max_index,i],:] #交换
        l = AugmentedMatrix[i+1:,i] / AugmentedMatrix[i,i] #计算系数
        m =  np.tile(AugmentedMatrix[i,:],(l.shape[0],1)) * np.tile(l,(col,1)).T #计算消元时减去的矩阵
        AugmentedMatrix[i+1:,:] = AugmentedMatrix[i+1:,:] - m #消元
    if(AugmentedMatrix[col - 2,col - 2] == 0):
            print("无唯一解")
            return
   #代入
    x = np.zeros(col-1)
    for i in range(col-2,-1,-1):
        x[i] = (AugmentedMatrix[i,-1] - np.dot(AugmentedMatrix[i,:-1] , x.T)) / AugmentedMatrix[i,i]
    return x    
#三角LU求解
def lu(scalex,scaley):
    m,n = scalex.shape
    L = np.zeros((n,n))
    U = np.zeros((m,m))
    #初始化L和U矩阵
    for i in range(0,n):
        if(scalex[i][i]==0):
            print("主元素为:")
    for i in range(0,len(scalex)):
        U[0][i] = scalex[0][i]
        L[i][0] = scalex[i][0] / U[0][0]
    #先计算U行，在计算L列
    for r in range(1,n):
        #i代表的是第r行的第i个参数
        for i in range(r,n):
            temp = 0
            for k in range(0,r):
                temp += L[r][k] * U[k][i]
            U[r][i] = scalex[r][i] - temp
            #计算第r行之后，计算第r列，但是r = n 时即主元素必为1
            #现在就是在r+1行就是i
            for i in range(r+1,n):
                temp2 = 0
                for k in range(0,r):
                  temp2 += L[i][k] * U[k][r]
                #print('temp2 = ',temp2)
                #print("L["+ str(i) + "]" + "[" + str(k) + "]" + " = ",L[i][k])
                #print("U[" + str(k) + "]" + "[" + str(r) + "]" + " = " ,U[k][r])
                #print(temp2)
                #print(L[i][k]," U" , U[k][i])
                #print(scalex[i][r])
                L[i][r] = (scalex[i][r] - temp2) /U[r][r]

               # print("L["+ str(i) + "]" + "[" + str(r) + "]" + " = " ,L[i][r])

        for i in range(1,n):
            if (r==i):
                L[r][i] = 1
    Y = np.zeros(n)
    X = np.zeros(n)
    try:
        Y[0] = scaley[0] / L[0][0]
    except:
        print(Y[0],scaley[0],L[0][0])
        print(type(Y[0]))
        print(type(scaley))
        print(type(L))
        exit()
    for i in range(1,n):
        for j in range(0,i):
            scaley[i]= scaley[i]- (Y[j]*L[i][j])
        Y[i] = scaley[i] / L[i][i]
    #print(Y)
    X[n-1] = Y[n-1] / U[n-1][n-1]
    for i in range(n-2,-1,-1):
        for j in range(i+1,n):
            Y[i]=Y[i]-(U[i][j]*X[j])
        X[i] = Y[i] / U[i][i]
    for i in range(0,n):
        print('x' +str(i)+ " = " ,X[i])
    #print(L)
    #print(U)
x = gaussinlist(scalex,scaley)
print(x)
scalex = np.array([[31.0,-13,0,0,0,-10,0,0,0],[-13,35,-9,0,-11,0,0,0,0],[0,-9,31,-10,0,0,0,0,0],[0,0,-10,79,-30,0,0,0,-9],
               [0,0,0,-30,57,-7,0,-5,0],[0,0,0,0,-7,47,-30,0,0,],[0,0,0,0,0,-30,41,0,0],[0,0,0,0,-5,0,0,27,-2],
              [0,0,0,-9,0,0,0,-2,29]])
scaley = np.array([[-15.0],[27],[-23],[0],[-20],[12],[-7],[7],[10]])
x = linalg.solve(radomx2,radomy2)
#print(x)
#testx = np.array([[1,1,7],[2,3,5],[4,2,6]])
#testy = np.array([2,3,4])

scalex2 = np.array([[30,33,-43,-11,-38,-29,37,28,23],[-480,-523,644,128,621,480,-618,-489,-329],
                   [60,266,-1862,-1991,464,546,-968,-1567,1652],[540,624,-782,290,-893,123,567,5,-122],
                    [-450,-675,2245,2326,-1512,1230,-822,129,-189],[-300,-120,-1114,-1295,1946,302,-376,-1540,-609],
                    [1080,998,508,2460,-1628,-1358,2896,2828,-2002],[-1080,-1408,3340,2267,21,-1202,866,-2690,-1351],
                    [-300,-435,1594,1685,340,2279,-27,2917,-2336]])
scaley2 = np.array([188,-3145,-4994,680,7845,1876,9712,-11599,10127])
#lu(radomx,radomy)












